#ifndef MATRIX_FREE_SOLVERS_HPP
#define MATRIX_FREE_SOLVERS_HPP
#include "HeaderFile.h"

//http://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html
namespace Chapter3_SparseLinearAlgebra
{

namespace Section3_MatrixFreeSolvers
{

// 诸如ConjugateGradient和BiCGSTAB之类的迭代求解器可以在无矩阵上下文中使用。
// 为此，用户必须提供一个继承EigenBase <>并实现以下方法的包装器类：

// Index rows()和Index cols()：分别返回行数和列数
// operator* 带有您的类型和Eigen密集列向量（其实际实现在internal :: generic_product_impl类的特殊化中进行）
// Eigen::internal::traits<> 还必须专门针对包装器类型。

// 这是包装Eigen :: SparseMatrix的完整示例：

// #include <iostream>
// #include <Eigen/Core>
// #include <Eigen/Dense>
// #include <Eigen/IterativeLinearSolvers>
// #include <unsupported/Eigen/IterativeSolvers>
// class MatrixReplacement;
// using Eigen::SparseMatrix;
// namespace Eigen {
// namespace internal {
//   // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
//   template<>
//   struct traits<MatrixReplacement> :  public Eigen::internal::traits<Eigen::SparseMatrix<double> >
//   {};
// }
// }
// // Example of a matrix-free wrapper from a user type to Eigen's compatible type
// // For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix.
// class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
// public:
//   // Required typedefs, constants, and method:
//   typedef double Scalar;
//   typedef double RealScalar;
//   typedef int StorageIndex;
//   enum {
//     ColsAtCompileTime = Eigen::Dynamic,
//     MaxColsAtCompileTime = Eigen::Dynamic,
//     IsRowMajor = false
//   };
//   Index rows() const { return mp_mat->rows(); }
//   Index cols() const { return mp_mat->cols(); }
//   template<typename Rhs>
//   Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
//     return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
//   }
//   // Custom API:
//   MatrixReplacement() : mp_mat(0) {}
//   void attachMyMatrix(const SparseMatrix<double> &mat) {
//     mp_mat = &mat;
//   }
//   const SparseMatrix<double> my_matrix() const { return *mp_mat; }
// private:
//   const SparseMatrix<double> *mp_mat;
// };
// // Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
// namespace Eigen {
// namespace internal {
//   template<typename Rhs>
//   struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
//   : generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
//   {
//     typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
//     template<typename Dest>
//     static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
//     {
//       // This method should implement "dst += alpha * lhs * rhs" inplace,
//       // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
//       assert(alpha==Scalar(1) && "scaling is not implemented");
//       EIGEN_ONLY_USED_FOR_DEBUG(alpha);
//       // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
//       // but let's do something fancier (and less efficient):
//       for(Index i=0; i<lhs.cols(); ++i)
//         dst += rhs(i) * lhs.my_matrix().col(i);
//     }
//   };
// }
// }
// int main()
// {
//   int n = 10;
//   Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
//   S = S.transpose()*S;
//   MatrixReplacement A;
//   A.attachMyMatrix(S);
//   Eigen::VectorXd b(n), x;
//   b.setRandom();
//   // Solve Ax = b using various iterative solver with matrix-free version:
//   {
//     Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
//     cg.compute(A);
//     x = cg.solve(b);
//     std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
//   }
//   {
//     Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg;
//     bicg.compute(A);
//     x = bicg.solve(b);
//     std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
//   }
//   {
//     Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
//     gmres.compute(A);
//     x = gmres.solve(b);
//     std::cout << "GMRES:    #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
//   }
//   {
//     Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
//     gmres.compute(A);
//     x = gmres.solve(b);
//     std::cout << "DGMRES:   #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
//   }
//   {
//     Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres;
//     minres.compute(A);
//     x = minres.solve(b);
//     std::cout << "MINRES:   #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
//   }
// }

// 输出：

// CG:       #iterations: 20, estimated error: 8.86333e-14
// BiCGSTAB: #iterations: 20, estimated error: 2.10809e-15
// GMRES:    #iterations: 10, estimated error: 0
// DGMRES:   #iterations: 20, estimated error: 1.10455e-28
// MINRES:   #iterations: 20, estimated error: 2.94473e-14
} // namespace Section3_MatrixFreeSolvers
} // namespace Chapter3_SparseLinearAlgebra
#endif